Squirrel

Download and load raster

The raster is freely available from the web

# Download globCover dataset
# download.file('http://due.esrin.esa.int/files/Globcover2009_V2.3_Global_.zip',destfile = "Globcover2009_V2.3_Global_.zip")
# unzip('Globcover2009_V2.3_Global_.zip',exdir = 'GlobCover')
# Load raster
glob = raster("GLOBCOVER_L4_200901_200912_V2.3.tif")

We also need to map the Globcover labels to a continous number, as indicated by the legend found here. We set water bodies to NA. Later on, we set the roughness ot water bodies to 100 to indicate that the squirrel cannot swim over them.

# map labels to roughness values
globCover = c(11,14,20,30,40,50,60,70,90,100,
              110,120,130,140,150,160,170,180,190,200,210,220,230)
roughness = c(0.10,  0.10,  0.30,  0.30,  1.50,  1.50,  1.50,  1.50,  1.50,  1.50,
              1.50,  0.50,  0.10,  0.03,  0.05,  0.50,  0.60,  0.20,  1.00,  NA,  NA,  NA,  NA)

Map of Spain

We chop a small section of the world map around the Iberian peninsula. The code is made to work only for squared rasters.

The start and end points of the path are set to be the northernmost point (Estaca do bares) and the sourthernmost one (Tarifa).

# Chop a small section
lat = 40
lon = -5.1
n=4.5 # width of the box
cropbox2 <-c(lon-n,lon+n,lat-n,lat+n) # crop a square box
spain_img = crop(glob,cropbox2)
# re-classify the raster values to numeric
small_img_plot = reclassify(spain_img,cbind(globCover,roughness))
# Let's tranform NA values to "-100" roughness so that we squirrel cannot swim. This one is only used for calculations.
small_img = reclassify(small_img_plot,c(NA,-100,-100))
# Save the rasters for later use in the dashboard
raster::writeRaster(spain_img, "Files/GlobCover_Spain.tif",overwrite=TRUE)
# Plot
# plot(spain_img) # Original globcover labels
plot(small_img_plot) # transformed labels, with NA
# plot(small_img) # we use this one for the analysis with NA=-100
# Start and end points
start_point = c(43.780824, -7.681479) # Estaca do bares
end_point = c(36.003797, -5.609124)# Tarifa
# Plot start and end points
points(x=start_point[2],y=start_point[1],pch=19)

points(x=end_point[2],y=end_point[1],pch=19)

Dimensions of the raster:

# translate raster to matrix
mat = raster::as.matrix(small_img)
# create 0-1 edge matrix
N = dim(mat)[1] * dim(mat)[2] # Number of vertices
Ni= dim(mat)[1] # number of rows
Nj= dim(mat)[2] # number of columns
N

[1] 10497600

Next we create the graph as a matrix. The matrix is equal to a 0 if two pixels are not linked, and equal to the roughness of theincoming pixel if they are. Two pixels are linked if they are adjacent to each other. This huge matrix of dimension 10497600 by 10497600 is really sparse and unless optimized takes too long to be created. It is build inside the function build_matrix_graph() and it took few tries to make it run at reasonable speeds.

The values of the raster are slighly modified to be used in a standard shortest-path algorithm: 0.0000001 weight on an edge means “forest”, while 1.5 means “bare land” and 100 means “water bodies”.

# Call our homemade function
val_raster = 1.5-values(small_img)+0.0000001
t = Sys.time()
mat_path = build_matrix_graph(Ni,Nj,val_raster = val_raster)
Sys.time() - t

Time difference of 32.85089 secs

mat_path[1:10,1:10]

10 x 10 sparse Matrix of class "dgCMatrix"
                                                                 
 [1,]   .   101.5   .     .     .     .     .     .     .     .  
 [2,] 101.5   .   101.5   .     .     .     .     .     .     .  
 [3,]   .   101.5   .   101.5   .     .     .     .     .     .  
 [4,]   .     .   101.5   .   101.5   .     .     .     .     .  
 [5,]   .     .     .   101.5   .   101.5   .     .     .     .  
 [6,]   .     .     .     .   101.5   .   101.5   .     .     .  
 [7,]   .     .     .     .     .   101.5   .   101.5   .     .  
 [8,]   .     .     .     .     .     .   101.5   .   101.5   .  
 [9,]   .     .     .     .     .     .     .   101.5   .   101.5
[10,]   .     .     .     .     .     .     .     .   101.5   .  

The shortest path is claulcated using the igraph library. It works surprisingly fast.

# Build igraph
g <- graph.adjacency(mat_path, mode="directed",weighted = T)
  
## Calculate shortest path between two poitns
# Find index of the start and end point
img_points = rasterToPoints(small_img)
idx_start = which.min( abs( img_points[,1] - start_point[2]) + abs( img_points[,2] - start_point[1]))
idx_end = which.min( abs( img_points[,1] - end_point[2]) + abs( img_points[,2] - end_point[1]))
                     
# shortest path algorithm
news.path <- shortest_paths(g,from = idx_start, 
                            to  = idx_end,
                            output = "both") # both path nodes and edges
# Find those poitns in lat,lon
path_squirrel = news.path$vpath[[1]] %>% as.numeric
path_points = img_points[path_squirrel,1:2]
# Save squirrel path
write.csv(data.frame(lon = path_points[,1],
           lat = path_points[,2],
           n = path_squirrel),
          file="Files/path_coordinates_solution.csv",row.names=F)

Visualize solution

The path and the roughness (land cover) map is displayed below. Due to its size, the land cover map needs to be aggregated by a factor of 3(i.e., reduce its size).